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Abstract 

The inclusive four-particle phase space integral of any 1^4 matrix element in massless QCD 
contains divergences due to the soft and collinear emission of up to two particles in the final state. We 
show that any term appearing in this phase space integral can be expressed as linear combination of 
only four master integrals. These four master integrals are all computed in dimensional regularisation 
up to their fourth order terms, relevant to next-to-next-to-leading order jet calculations, both in an 
analytic form and purely numerically. New analytical and numerical techniques are developed in this 
context. We introduce the tripole parametrisation of the four-parton phase space. Furthermore, we 
exploit unitarity relations between multi-parton phase space integrals and multi-loop integrals. For 
the numerical calculation, the iterated sector decomposition of loop integrals is extended to phase 
space integrals. The results in this paper lead to infrared subtraction terms needed for the double real 
radiation contributions to jet physics in e^e~ annihilation at the next-to-next-to-leading order. 



1 Introduction 



In recent years, data on jet observables from LEP, HERA and the Tcvatron have reached an impressive 
level of experimental precision. Using these data for a determination of the strong coupling constant, or 
as a measurement of parton distributions in the colliding hadrons, it turns out that the error on these 
extractions is largely dominated by the uncertainty on the underlying theoretical calculations, which are 
accurate only to the next-to-leading order (NLO) in perturbative QCD. 

In order to improve on this situation, theoretical calculations of jet observables accurate to the next- 
to-next-to-leading order (NNLO) are mandatory. For each observable, they require several ingredients: To 
compute the corrections to an n-jet observable, one needs the two-loop n-parton matrix elements, the one- 
loop (n-l-l)-parton matrix elements and the tree level (n-|-2)-parton matrix elements. In the recent past, 
enormous progress has been made especially on the calculation of two-loop 2 — > 2 and 1 — > 3 QCD matrix 
elements, which are now known for all massless parton-parton scattering processes ^EUHlll] relevant to 
hadron colliders as well as for 7* qqg and its crossings [7] . Most of these results were even verified 
independently by different groups. For the corresponding partonic processes, the one-loop matrix elements 
with one additional parton and the tree-level matrix elements with two more partons are also known and 
form part of NLO programs for 1 — > 4 |S] and 2^3 reactions 0. Since these matrix elements lead 
to infrared singularities due to one or two partons becoming theoretically unresolved (soft or collinear), 
one needs to find one- and two-particle subtraction terms, which account for these singularities in the 
matrix elements, and are sufficiently simple to be integrated analytically over the unresolved phase space. 
One-particle subtraction at tree level is well understood from NLO calculations ^1^2^] and general 
algorithms are available for one-particle subtraction at one loop |13II14[[T5] . in a form that could recently 
be integrated analytically [T^IT^ . Tree- level two-particle subtraction terms have been extensively studied 
in the literature jl()lll7l[THj . their integration over the unresolved phase space could however up to now 
be performed only in one particular infrared subtraction scheme (the hybrid subtraction method) in the 
calculation of higher order corrections to the photon-plus-one-jet rate in e+e~ annihilation |19[l2l)j . The 
same techniques (and the same scheme) were used very recently in the rederivation of the time-like gluon- 
to-gluon splitting function from splitting amplitudes A general two-particle subtraction procedure, 

consisting of subtraction terms which can be safely implemented numerically and at the same time are 
sufficiently simple to be integrated analytically is still lacking at the moment. It is the aim of the present 
paper to contribute to such a method by classifying and computing all four-particle phase space integrals 
relevant to 1 ^ 4 parton reactions in massless QCD. 

One of the most widely used infrared subtraction schemes at NLO is the dipole formalism established 
by Catani and Seymour |12j . The formalism relies on an exact factorisation of the n-particle phase space 
into an (n~l)-particle phase space and a so-called dipole phase space, which is, up to a normalisation factor, 
equal to a massless three-particle phase space. In this formalism, one derives (process-independent) dipole 
subtraction terms for all single unresolved partonic configurations, which are then integrated over the dipole 
phase space (corresponding to a three-parton configuration of emitter, unresolved parton and spectator). If 
combined appropriately, these integrals of the subtraction terms over the dipole phase space can be related 
to integrals of 1 3 tree level QCD matrix elements integrated over the three-particle phase space. 

At NNLO, one encounters configurations of two unresolved partons between emitter and spectator. 
As will be shown below in Section 14.11 these can be accommodated into a tripole formulation, relying 
on the exact factorisation of the n-particle phase space into an (n-2)-particle phase space and a tripole 
phase space. This tripole phase space is proportional to the four-particle phase space. Based on this 
observation, one might envisage that all double unresolved partonic configurations which correspond to a 
single emitter-spectator pair may be expressed by appropriate tripole subtraction terms. To facilitate their 
analytic integration over the tripole phase space, these may be further combined to relate to integrals of 
1^4 tree-level QCD matrix elements integrated over the four-particle phase space. In this paper, we 
show that all integrals of this type can be expressed as linear combinations of only four so-called master 
integrals, which we compute both analytically and numerically. 

The reduction to master integrals follows largely methods developed for multi-loop integrals, and is 
briefly presented in Section |21 To compute the resulting master integrals analytically in Section 01 we 
develop first the tripole formulation of the four-particle phase space, which is then used to calculate two of 



1 



the master integrals. We then show that the remaining two master integrals can be inferred from unitarity 
relations between known multi-loop integrals and phase space integrals. Further, one of the master integrals 
calculated explicitly, as well as a reducible integral, are rederived from unitarity relations as crosschecks. 

In Section |S1 the master integrals are calculated numerically. The infrared poles are isolated by an 
automated sector decomposition algorithm which already has proven useful in the calculation of multi-loop 
integrals |22II23| . and which has been extended in 24 and in the present work to tackle infrared divergent 
integrals over a multi-parton phase space as well. 

A summary and conclusions are given in Section ^ The Appendix contains formulae relevant to the 
massless (1 — > n)-parton phase space. 



2 Kinematics and Notation 

In this paper, we consider the decay of a massive particle Q (Z-boson, off-shell photon or Higgs-boson) into 
four massless QCD partons (Pi, . . . P4: qqgg, qqqq or gggg): 

Q{q) ^ PiiPi) + P2{P2) + P-siPs) + PiiPi) , (2.1) 

where q,pi, . . . ,P4 denote the momenta of the particles. The kinematics of this process are uniquely defined 
by the six Lorentz invariants Sij = 2pi ■ pj formed by any pairs of partons. These appear in the numerator 
and denominator (propagators) of QCD matrix elements. Only five of them are independent, the sixth 
being fixed by energy-momentum conservation 

q^ = S12 + Si3 + Si4 -I- S23 + S24 + S34 • 
Besides these, QCD matrix elements also contain propagators yielding triple invariants 

Sijk — Sij S^/c 4^ Sjk • 

Without loss of generality, we can always choose Pi and P2 such that the associated scalar product S12 and 
triple invariants S123, S124 do not appear as a propagator in a given term of the squared matrix element. 
Any single term in a squared four-particle QCD matrix element can then be expressed as 

T{S-^ . . . S-^;D'r . ■ ■ DD = j^Lr^^rn, (2.2) 

where the Si are the Lorentz invariants Sij and the Di are the propagators, being either invariants Sij or 
triple invariants Syfe. The rii and mi are positive integers. It is assumed that whenever possible, scalar 
products in the numerator are expressed as linear combinations of the propagators, such that only a minimal 
number of different scalar products is present in the numerator. In fact, one finds that for t — 4 different 
propagators, only 9 — t different irreducible scalar products exist. The integral of the term 1)2. 2|l over the 
four-particle phase space is denoted by 

It (ST . . . S;'" ; D"^'- . . . D'^' ) = J dPS^ TiS"^' . . . S^'" ; . . . D'^' ) . (2.3) 

The phase space measure dPS'4 is defined in (|A.5p in the appendix. 

The topology of It (interconnection of propagators and external momenta) is uniquely determined by 
the set of propagators {D^, . . . ,Dt). In order to visualise the topology of It in the form of a cut diagram, 
we use the Cutkosky rules |25ll2f)j to introduce the four cut-propagators 

^-2...-(,D = ^-^ (-1,...4). (2.4) 
With the help of the cut-propagators, the four-particle phase space measure in (I2.3|) can be written as 

id -id -id id -1 

mo .-4 dpi d P2 dps d P4 sd cdf \ 1 rn t:\ 

^^^^ = * (2^ (2^ (2^ (2^ ^(<Z-Pi-I>2~P3-P4) ^^^^^3^^ ^ (2-5) 
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such that the integral It can be expressed as a cut through a three-loop vacuum polarisation diagram. Any 
integral It for a given topology (fixed set of i — 4 propagators) can be attributed to a class of integrals 
lT{t,r,s) of identical numerator mass dimension s = X]i=i ^^'^ denominator mass dimension r = Yli=i 
It should be noted that we also allow the cut-propagators to appear in general with integer powers larger 
than unity, without further specifying the physical meaning of those. From combinatorics, one finds that 
the class lT(t,r.s) contains 

'r - 1\ /8 - i + 

different integrals. 



Nt,r,s 



3 Reduction to Master Integrals 

The number N{It.r,s) of the integrals grows quickly as r,s increase, but the integrals are not all linearly 
independent. These linear relations are the so-called integration-by-part (IBP) identities, which were orig- 
inally formulated for multi-loop two-point functions |27II28| with ordinary loop propagators, but can be 
safely extended to multi-loop integrals with cut propagators j^H- These identities follow from the fact 
that the integral over the total derivative with respect to any loop momentum vanishes in dimensional 
regularisation 

where J is any combination of propagators, scalar products and loop momentum vectors. J can be a 
vector or tensor of any rank. The IBP identities for three-loop two-point functions were derived and 
solved symbolically a long time ago |2H| in the context of the calculation of the three-loop /3-function. The 
solution of these identities was implemented into the program package Mincer |29| written in the algebraic 
programming language Form |2Q||. Mincer was widely used in recent years to compute a large number of 
multi-loop corrections which correspond to three-loop two-point functions or can be related to them by 
expansion. Reviews of these results can be found in |31[l32j . 

In recent times, IBP-type relations between multi-loop integrals with more than two legs have been 
studied extensively. In particular, it was found that scalar multi-loop integrals with three or more legs 
fulfil, besides the ordinary IBP relations, also relations following from the invariance of the integrals under 
an infinitesimal Lorentz transformation. These Lorentz invariance (LI) identities, combined with the IBP 
identities, were used extensively in the computation of two-loop four-point matrix elements ^ElElEElini- 
In a related development, it was observed that the IBP identities for Feynman integrals with the same total 
number of external and loop momenta are equivalent 32, , a development that was used soon thereafter to 
compute two-loop three-point functions j34| . 

Exploiting the observation that differentiations involved in the integration-by-parts identities are insen- 
sitive to the imaginary parts of the two terms in a cut-propagator 122], one can generalise the IBP method 
to phase space integrals over tree-level or loop matrix elements. In the solution of the IBP equations for a 
given cut integral, one obtains integrals with one or more of the cut propagators eliminated. These inte- 
grals can be discarded, since they do not contribute anymore to the same cut of the underlying multi-loop 
diagram. 

To derive the reduction of any given integral It of the form (|2.3|l . we use the same strategy |^ as 
applied in the reduction of the two-loop four-point functions in ,36, by generating all IBP identities for 
integrals up to the maximum values of r and s required for the integration of QCD four-particle matrix 
elements. From dimensional counting, one finds r < 8, s < 3. The solution algorithm j35| proceeds by 
solving all identities up to these values using the computer algebra languages Form and Maple ^Tj . 

After carrying out the reduction for all topologies that can be formed from the invariants and triple 
invariants involving the four final state momenta, we find that all inclusive four-particle phase space integrals 
of massless QCD matrix elements can be expressed as a linear combination of four master integrals. These 
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four master integrals are: 




The reduction procedure generates coefficients containing explicit poles of up to 1/ in front of R4 and up 
to in front of Rq, while the coefficients multiplying Rs,a and i?8,& are finite. 

In the following two sections, we will compute these four master integrals both analytically and numer- 
ically to a sufficiently high order in e required for the calculation of integrated four-particle QCD matrix 
elements. 



4 Analytic Calculation of Master Integrals 
4.1 Tripole Phase Space 

The dipole formalism |12j of Catani and Seymour is a subtraction formalism for the calculation of jet 
cross sections at NLO. The main features of this formalism are factorisation formulae for matrix elements 
and phase space which locally reproduce single soft and single collinear behaviour of the cross sections 
in the appropriate limits. The corresponding dipole contributions to the cross sections can be integrated 
analytically and numerically over the whole phase space. In particular, the exact factorisation of the three- 
parton phase space into a dipole phase space and a two-parton phase space is obtained by redefining a 
set of three massless on-shell momenta (emitter, unresolved parton, spectator) into two on-shell massless 
momenta. This procedure was further generalised by Kosower |17| , also accounting for the recombination of 
more than three partonic momenta. In |17| a non-linear remapping of momenta, which smoothly interpolates 
between all potentially singular regions (as required for the numerical implementation) was introduced. 
Such a non-linear remapping is however inappropriate for analytic integration, and we shall derive a linear 
remapping below. 

As an extension of jl2| . we introduce a parametrisation of momenta allowing the factorisation of a 
massless n-parton phase space into a massless (n-2)-parton phase space times a tripole phase space factor: 

dPSnipi, ■ • ■ ,P«) = dPSn-2iPl, ■ • ■ ,Pn-4,Pn~3iPn-2)dPST , (4.1) 

In particular, it enables to factorise exactly the phase space of four on-shell parton momenta (pi..p4), 
which can also denote the four-particle contribution to an (n > 4)-particle phase space, into a two-parton 
phase space times a tripole phase space factor 

dP54(pi, . . . ,P4) = dPS2{p~^,P2) dPSr. (4.2) 

In this parametrisation, pi is the momentum of the emitter, p2 the one of the spectator, while p^ and p4 
are the momenta of the unresolved partons. 

The tripole momenta {p^^,P2 ) are combinations of four on-shell parton momenta. They are defined 

by 

Li LL.a.a 2/134,2 u 

= ^ — p^ (4.3) 

1 - yi34,2 
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where the dimensionless variable j/134,2 is 

P1-P3 + P1-P4 + P3-P4 f. .^ 

2/134 2 — • (4.4) 

PI.PZ + Pl.pA + P3 -Pi+Pl ■P2 + P2 -Pa + P2 -Pi 

The tripole momenta are on-shell (p~ — p\ — ^ ) and momentum conservation is implemented exactly 

Pi + P2 + P3 + K = P~ + P2 ■ 
To derive the tripole factorisation of dPS'4, we start from its definition in terms of the invariants 

dPSi = (2^)4"3'i (q2)3-f 2l-2<i(_A^)^ 0(„A4) ^(si2 + S13 + 514 + S23 + S24 + 534 - 9') 
AVld-l drJd_3dsi2 dsi3 dsi4 dS23 dS24 ds34 , 

where the four-particle Gram determinant A4 is given by 

I A 2 2 I 2 2 I 2 2 

+ A4 = S12 S34 + Sl3 S24 + Sl4 S23 



-2|^Si2S23S34Sl4 + S13S23S24S14 + S12S24S34S13 
= A (S12S34, S13S24, S14S23) (4.5) 

and A is the Kallen function, A(a:;, y, z) — -\- + — 2{xy + xz + yz). 

In order to obtain the massless four-parton phase space di-'S'4 in terms of dimensionless invariants 
defined with P^iP2 ^ given in 14.21) . the invariants Sij need to be substituted as follows 

523 = (1 - 2/134,2)523 = (1 ~ 2/134) ^1 , 

524 = (1 - 2/l34,2)s24 = (1 - 2/134) Z2 , 

512 = (1 - 2/l34,2)Si2 = (1 - 2/134) (1 - ^1 - Z2) . (4.6) 

The last simplification occurs since, in the case of only four partons in the final state, one has 

^P^rP2 = ^{^4 = ^1234 = (4.7) 

and j/134,2 is the dimensionless variable j/134 — 5134/9^. The fractions of momenta Zi are defined by: 

P3-P2 



Zl 



PlT4-P2 



z, = J^. (4.8) 

In the tripole phase space dPSr, the variables S13, S14 and S134 remain unchanged and the triple invariant 
S134 is used instead of S34. 

From the four-parton phase space, one can easily factorise a two-parton phase space factor dP52 as 
defined in the Appendix using (|4.7I) . The integration variables in dPS'4 become 

dsi2 dsi3 dsi4 ds23 ds24 dsi34 = (1 - 2/134)^ (q'^f dsi234 dyi^ dyi4 dzi dzi d2/i34 , 

where the dsi234 on the right-hand side becomes part of the two-parton phase space while the other inte- 
gration variables contribute to the tripole phase space factor dPSx- Furthermore, the Gram determinant 
A4 is rescaled as 

A4 - (1 - 2/134)' {q'r a; , 

with 

A4 = A (2/12(2/134 - 2/13 - yw), yi3Z2, yuzi) . (4.9) 
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The factorised form of the four-parton phase space yielding the tripole phase space factor then reads: 

dPSi = dPS2 dPSr , (4.10) 

with 

dPSr = {q' f-'' r(l-2e) ~ 2/i34)^-'^(-A^)-i/2- e(-A;) dyi3 dyu dz, dz^ dym ■ (4.11) 

This tripole phase space factor dPSx is proportional to the four-parton phase space, since the two- 
parton phase space is a constant. It is also proportional to the triple coUinear phase space factor given 
in [201 obtained by requiring a special kinematical situation in which the three partons whose momenta are 
Pi,P3,P4 are coUinear. In [211 this triple coUinear phase space was used to compute the integral over the 
subtraction term appropriate to this kinematics, the triple coUinear splitting function. The normalisation 
factor relating triple coUinear phase space and tripole phase space is just (1 — 2/134)^"^"^. It follows clearly 
that in the triple coUinear limit where j/134 0, the tripole phase space factorisation is appropriate. 

The volume of the four-parton phase space is well-known. However, we shall now rederive it using the 
tripole parametrisation described above, to demonstrate an application of the factorisation formula H4.10|) 
and to illustrate the necessary variable transformations. 

In this formula, — A4 can be written as a quadratic in j/13, the first variable to be integrated over. We 
obtain 

-A4 = (1 - Zif{yi3 - yi3,a){yi3,b - yi3) , 
with 2/13, a, 2/13,6 being the roots of the quadratic. Using 

2/13 = (2/13, f) - 2/13, a)x + 2/l3,a, 

one has 

I (-A;)-i/2- d2/i3 = (1 - z,)-'-'^ iy,3,b - 2/i3,a)-'^ £ dx [x(l - xr'^'-' ■ 
The root difference (2/13.6 — 2/13. a) satisfies the relation 

(1 - Zi)" (2/13, b - 2/13, a)^ = 162:1 Z22/14 {l~ Zi~ Z2) (2/134(1 - Zl) - 2/14) = <^yi3- 

Requiring that the Gram determinant — A4 is positive suggests the following change of variables 

2/14 = 2/134(1 - zi)v , 

Z2 = (1 - Zl)/ . 

with t,v ranging between and 1. Then Sy-^^ completely factorises in these new variables: 

Syi3 = 16 2/m (1 - ^i)" zit{l-t)v{l-v) , 
such that the volume of the four-parton phase space then becomes: 

dPS'4 = / dPS2 {q^f-^' 2-1" 7r-5+2^ ^ 



r(l-2e) 



^'(1/2 e) ^^^^^^^ _ ^^^^^,^2. r (1 _ ^^)l-2e 



r(i-26) Jo 







At t-' (1 - ty [ dv (1 - v)-'. (4.12) 

JO 

The integrations can be performed straightforwardly. Finally the volume of the four-parton phase space 
rScLcis ' 

7? - P - /dPS - P . 2.2-2e r3(l-,)r(2-26) 

R^-P^-J ^PS,-P2 ( q ) r(3-3e)r(4-4e) ' 
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where the volume of the two-particle phase space P2 lA.9p has been factored out. Introducing the common 
normalisation factor 



Sv = P2 



167r2r(l-e) 



one arrives at 
i?4 



2 r5(l-6)r(2-2e) 
' r(3-3e)r(4-4e) 



1 59 



/2263 



/ 72023 597r2 13(3 



12 72 I 432 9 / ' V 2592 54 6 



/ 2073631 22637r2 767C3 



V 15552 324 36 1080 



(4.13) 



(4.14) 



(4.15) 



4.2 Calculation of R 



Using the factorisation of dPS4 into the tripole phase space dPSx and P2 — J dPS2, the master integral 
i?8,a defined by 

i?8,a- fdPSi- ^ 



can be written as follows: 



R8,a = P2 J dPSril-yisirHq 



2 i'„2\-i 



Sl3 Si4 S23 S24 



1 

yi3 yi4 zi Z2 



P2 {q^)-^-^' 2-1° 7r-5+2^ 



r(l-2e) 

/(I - 2/i34)-'-'^("A;)-i/2-^e(_A;) dyi3 dyi4 dzi dz2 dyi34 i 

J yi3 yu 



yi3 yu zi Z2 



(4.16) 



To perform the integrations, one factorises the Gram determinant using the same change of variables 
as for the evaluation of R4. However due to the presence of the invariants in the denominator of the 
integrands, the integrations are more complicated. Indeed, hypergeometric functions of involved rational 
arguments occur at different stages of these five subsequent analytic integrations. Using several non-linear 
transformations |38ll39j on these functions where appropriate, all integrals can be carried out in a closed 
analytic form, yielding finally 



Rs,a = Sr [q 



6 r5(l-e)r(l-2e) 



£4 r(l-3e)r(l-4e) 



1^3 (1, -e, -e, -e; -3e, 1 - e, 1 + e; 1) 



1 r3(l-e)r(l + e)r3(l-2e) 



r2(l-4e) 



;i^2(-2e,-2e,-2e;-4e,l-2e;l) 



Expanded up to finite order in e, i?8,a reads 



Rs.a = Sr {q ) 



2\-2-2e 



207r2 



I26C3 7^ 
e 18 



Oie) 



(4.17) 



(4.18) 



4.3 Unitarity Relation and Calculation of Rq and Rs^ 



The use of unitarity to obtain relations between loop integrals and phase space integrals can be found 
in various contexts in the literature. In |40) . unitarity is exploited to compute one- loop integrals from 
known single unresolved phase space integrals. To do the inverse, i.e. to derive phase space integrals 
from multi-loop integrals, has been exploited for example in |41) . In j42j. it was moreover attempted to 
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implement the unitarity cancellation of infrared singularities to compute jet observables purely numerically 
as appropriately weighted loop integrals. 

As outlined in Section|21 one can view the phase space integrals as particular cuts of multi-loop diagrams. 
In the case considered here, the four-particle phase space master integrals are particular cuts of massless 
three-loop two-point functions, which were studied extensively in the literature |43II28[I^ . In particular, 
all master integrals appearing in these three-loop two-point functions are either known exactly or to a very 
high order in e. These three master integrals are 43..28) : 




(4^)3' r(l + 3e)r(l-36)r4(l-e) ^ ^ 



(16^2)3 3er(3-3e)r(4-4e) 



(4.19) 



/r = 




(47r)3' r^(l-e)r^(l + e) 



(16^2)3 3e3r3(2-26) 



l + e + e2 + e3(i4(3_7) 



-67 + WCs 



2l7r^ 

"go" 



0(e5 



2\-3e 




0(e°) 



(4.20) 



(4.21) 




The optical theorem relates the imaginary part of a loop diagram to the sum over all its cuts in the 
form of a unitarity relation: 

2 Im = 

where i enumerates all possible ways to cut the vacuum polarisation diagram into two connected amplitudes 
under the condition that the final states resulting from the cut lines belong to physical processes. It has to 
be kept in mind that this equation relates matrix elements, not just loop and phase space integrals. In the 
case of the master integrals considered here, one obtains the corresponding matrix elements by considering 
a scalar theory involving scalar three- and four-point vertices and two-point as well as four-point couplings 
to an external scalar current. Accordingly, one has to include a factor i for each propagator and a factor 
—i for each vertex to obtain the desired relations among master integrals. 

If all but one of the cuts are known, the unitarity relation can be used to infer the remaining unknown 
cut. In the case considered here, the three-loop two-point functions can have 

1. two-particle cuts: two- loop vertex function or product of two one-loop vertex functions integrated 
over two-particle phase space. 

2. three-particle cuts: one-loop four-point function integrated over three-particle phase space. 

3. four-particle cuts: tree-level matrix element integrated over four-parton phase space. 

As we shall see below, the functions appearing in the two-particle cuts are well known from other calcula- 
tions, and the three-particle cuts can be computed more easily than the four-particle cuts. The unitarity 
relation following from the optical theorem can therefore be used to determine the remaining four-particle 
phase space master integrals Rq and Rg^b- Before computing these, we illustrate the application of the 
method on the rederivation of i?4. 

The unitarity relation relevant to i?4 reads pictorially: 



2 Im 




which yields 



2Im/4 = i?4. 



(4.22) 



Im ( (-1)""' r(l + ne) T{1 - ne) — ] ^ n 



Using the relation 

which is valid to all orders in e, one can extract R4 from the value of I4 listed in (|4.19|l : 



(4.23) 



r(3-3e)r(4-4e) 



in accordance with the result obtained by explicit computation in H4.14|l above. 
To determine Rg, we use the unitarity relation obtained for /g: 



2 Im 




reading 



= 2 Re 



2 Im /« = -2 Po RcAa +2 Re 



(4.24) 



The two- loop vertex integral A4 is known from form-factor calculations |44j . and appeared also in the 
reduction of on-shell and off-shell massless two-loop four-point functions |45j : 



A4 = 



167r2 



r(l - 26) 



1 5 /19 

2? ^ 2^ ^ vT ^ y 

211 197r2 
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65 Stt" 
IOC3 ) e^+0{e^) 



- 2C3 U 



Therefore, Rq can be immediately read off from (|4.24|) to be 
i?6 = Sriq^)-^' 



-1 + Y + M-12+ — + 9C3 



0(e^ 



(4.25) 



In order to obtain R^^b, we employ 



2 Im 




= 2 Re 




4 Re 




2 Im /g = -2 P2 Re As -I- 4 Re Vg 4- i?8,a + 4i?8,6 ■ (4.26) 
Again, the crossed vertex master integral is well known from the literature |44[I45| to be 



A, = 



{AttY r(l + e)r2(l-e) 



167r2 



r(l-2e) 



-2-2e 



1 57r-^ 23^ 1037r4 ^ 
6e^ e 180 



The three-particle cut Vs is the integral of the one-loop four-point functions over the three-particle 
phase space 

Vs = -i [ dPSs — Box(si3, S13, S12) , (4.27) 

J Si2 

with the one-loop box taking the well-known 44 compact form 



Box(si3,S23,Si2) = 



{AnY r(l + e)r2(l-e) 



167r2 



r(l - 2e) 



-2-e 2i 1 



{^1 > 72 



9 



yi3y23 
1 - yrs 

yi32/23 



1 



2/23. 
V13V23, 



2F1 \^-e,-e; 1 - e; 
2F1 i -e, -e; 1 - e; 



(1 - ?/i3)(l - ^23) 



1 - yi3 - 2/23 

1 - yi3 

1 - yi3 - ;/23 
1 - y23 



1 - yi3 - 2/23 



(1 - 2/i3)(l - ^23) 



where in this case = S12 + S13 + S23 and ?/y = Sijjq^. To perform the integration, the three-particle 
phase space is most conveniently parametrised in terms of the yij (see Appendix): 



dPS^ = (27r) 



-5+4e r,-5+2£ /„2Nl-2e 



- dri3_2edri2-2e (2/122/132/23) ' dyi2 d^is dy23 (1 - 2/12 - 2/13 - 2/23) 



Similar integrals were considered in the literature in the context of the integration of subtraction terms 
appropriate to single unresolved limits of one-loop amplitudes [THIT3] . While the first two terms can in fact 
be straightforwardly integrated over the three-particle phase space, more care is required for the last term. 
We have used several analytic continuations of hypergeometric functions |38| at intermediate stages where 
appropriate. Furthermore, some of the integrals could only be carried out on the series representation of 
hypergeometric functions, the resulting expressions were resummed using Mathematica '46'. One finally 
arrives at 

5 , 89^^ 13^4 



Re = Sr [q^y^-^' 
As a result, one infers from H4.26|l : 

-^8,6 — Sriq^y'^^'^'^ 



2e4 2e2 



180 



3 

4? 



177r2 44C3 



12e2 



eiTT^ 



-0{e) 



(4.28) 



(4.29) 



Finally, we can also employ the unitarity relation to check the correctness of the reduction developed in 
Section |31 by computing a reducible integral directly through the unitarity relation and comparing it with 
the result obtained from the reduction procedure. For this exercise, we choose 



i?8 




= / dPSi 



1 



S135134S24S234 



(4.30) 



which reduces to 



R. 



•8,r 



(1 _ 76+ 18e2) (2 - 3 6) (3-4e) 3 (1 - 26) (1 - 36) 

-fl4 H —5 Kg 



= Sv{q^r 



26^ 

1 

4? 



1262 



7C3 
6 



237r4 
^5" 



262 

-0{e) 



This integral is related to the imaginary part of |43[l^ 



h,r — ■ 




.= 0(6°), 



(4.31) 



(4.32) 



for which the unitarity relation reads 



2 Im 




= 2 Re 



4Rc- 




2Im/8,r = -2P2Re A6,..-h4Rey8,r- + |v43,r|^-P2 + 2i?; 



•8,r 



(4.33) 
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The terms in this relation are easily evaluated as 



A. 



167r2 

2\-2-2e 



(4^)^ r(l + e)r2(l~e) 
r(l-2e) 

2e4 + 6e2 



5C3 



5^ 
36 



J__^_ 6 
4e4 4e2 g''^ 



0(6) 



20 



I A 



3,r| 



(4^)^ r(l + e)r2(l-e^ 
167r2 



r(l-2e) 



/ 2\-2-2e 1 
) ^ 



such that together with H4.31|l . the unitarity relation H4.3r{|) is fulfilled identically. Since only the evalua- 
tion of i?8,r relied on the reduction of phase space integrals, this provides a strong check on the correct 
implementation of the reduction algorithm described in Section |2| 



5 Numerical Calculation of Master Integrals 

As outlined in 24 , the method of sector decomposition 47 to isolate the (infrared) divergences from a 
given parameter integral and to calculate the resulting pole coefficients numerically can be applied not 
only to multi-loop integrals j22ll2Hj . but also to phase space integrals. To this aim, the phase space for the 
1 — > 4 particle reaction considered here has to be cast into a convenient parameter integral form, as will be 
explained below. 

The starting point is the expression HA.5() for the four-particle phase space. This form, where the 
integration variables are the kinematic invariants s^, is most convenient for the sector decomposition 
algorithm since it is maximally symmetric in the integration variables. In addition, the potential infrared 
divergences stem a priori only from the limits Sij ^ of a certain set of invariants (appearing in the 
denominator of some matrix element), that is, they are located at the origin of parameter space. 

Further, we rescale the Mandelstam invariants, in order to deal with dimensionless parameters, by the 
definitions 



to arrive at 



6 



I dPSi = Cr(<z')*-^ /"{nda;, e(^,)}<5(l~^a;,0 [-A(siX6,X2a;5,X3a:4)] ' e(-A) , (5.1) 

j=i i=i 

where 

Cr = {2TT)^-^'^V{d - l)Vid - 2)V{d - 3) 2^~^'^ , 

A(a;ia;6, x^x^) is the Kallen function already introduced in (|4.5|) and V{d) is the c?-dimensional volume 
defined in HA.2|) . 

After these transformations, we obtain for the integrals Rq and Rsa,b 



Re = {q')-' APSi- -, (5.2) 

J (X2 + X4 + X6)[X3 + X5 + Xe ) 
R8,a = {qY^ fdPSi , (5.3) 

J X2X3X4X5 



R8,b = (g2)-4/dP54 . ^ ^ V ^ a- ^ ■ (5-4) 

J X2X3{X2+X4+Xe)[X3+X5+Xe) 

6 

The variables in l|5.1|l are constrained not only by the momentum conserving delta-distribution (5(1— ^ Xi), 

i=l 

but also by the requirement 

X{xiXq, X2X5, a;3a;4) < 4^ x^Xq + x\x\ -f x\x\ — 2 [xiX3X4^X(, + X2X3X4^x^ + xiX2X5Xq) < . (5.5) 
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Therefore the integration variables are hard to factorise and disentangling the overlapping pole structure 
in general is a highly nontrivial task. That is where the iterated decomposition into sectors until a non- 
overlapping form is reached shows its virtues. 



5.1 Sector decomposition 

First we decompose the original integral into a sum of six integrals corresponding to the six "primary 
sectors" created by inserting unity in the form 

6 6 

1 = E n - • 

j = l k=l 

In each primary sector j created in this way we transform variables according to 

_ ( Xj tk lik^j 
' I if fc = J 

and integrate out Xj using the momentum conserving ^-distribution, to arrive at the following form for the 
phase space integral (where we choose j = 1 to give a concrete example) : 

/ dPS[-'^ = Cr (g^)--^{ n / Q{U) 6(1 - i.)} 4"/.^) 

d-5 

[- \^'"'^\h,t2hMtA)\ ' e(-A(^^^i)), (5.6) 



A(sec\) = 1 + 



6 



i=2 



-X^'^^^XhMhMU) = tl + tltl + tltl-2{t3UtQ + t2hUU + t2hh) ■ (5.7) 

Solving the equation — A'^'^'^^^ = for tg (or, in general, in sector j, for ty-j, which occurs simply as t^-j 
in A("'='=J)): 

= ^2^5 + ^3^4 ± 2\/i2i3^4^5 — (^^2^5 ± V^3^4)^ 



and substituting 

where now 
leads to 



te — ie {tt ^ ^6 ) + ^6 ^ ^6 ^"sfh^t^t^ + ^2^5 + ^3^4 — '2^/t2t^t4tl , (5.8) 



d?6 [^6(1 - h)]--^^' e(t6) 9(1 - [i, (4 - tg-) + 1^]) . (5.9) 



^ 2 ^' i=2 ■ 

^-4+4c 
(seel) 

However, as a consequence of the transformation H5.8|l . singularities at = 1 can also occur now. For exam- 
ple, a factor l/tg that might be present in the matrix element can become singular for Iq — > 0, ti=2,3,4,5 1- 
As the subtractions performed later rely on the fact that only the limits ti Q can lead to a singularity, these 
potential singularities are remapped to the origin in the following way: The program checks if a (non-empty) 

set T exists such that a denominator vanishes if all ti G T are one and 0. For the U € T, it is checked if 

1-1 

in addition ^ is singular. If yes, the integration region is split according to dti — dti + Ji dti and 
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the resulting integrals are remapped to integrals from zero to one. Otherwise, simply the transformations 
ti ^ 1 ~ ti are performed for the U e T. Note that the limits {ik 0,ti ^ 1 [ti e T)} do in most cases 
not lead to a real 1/e singularity upon integration, even if they make the denominator vanish, but they will 
be remapped nevertheless for the sake of numerical stability. 

These transformations as well as the sector decomposition are iterated if necessary, that is, if a (com- 
posite"'^) denominator of the integrand still can vanish at a certain phase space point. 

The iterated sector decomposition proceeds as follows: The program determines the minimal set S ~ 
{tai , . ■ . , ^Qr} of parameters which in the limit ict — * make a composite denominator of the matrix element 
vanish. (The parameters t have been renamed t). Then the corresponding integrals are decomposed into r 
subsectors according to 

r r r 

n e(i - i„j = 5] nQ(^"'= (5.10) 

i—l k=l i=l 

In each new subsector, the variables are remapped to the unit cube by the substitution 

, _^ / taja^ for i^k , , 

\ ia, for z = • ^^-^^^ 

By construction, ta^, factorises from at least one of the denominator functions. The sector decomposition 
procedure is iterated until the denominators are of the form D = 1 + f{ta), where f{ta) is a non-negative 
function of the parameters ta- Then subtractions of the singularities, which are now all factorised explicitly, 
are carried out, using identities like 

/ <itata^~''^'^ ^{ta,tf3^a) — •?^(0, t/J^^a) + / dta t~^~ '^'^ I T {ta, tp^a) ~ ^{0, t p^a) 

JO ^£ Jo I 

where liint^^Q J-'(ta,tfj^a) is finite by construction. 

After having carried out the subtractions for each ta, the expansion in e can be performed. The result 
is a Laurent series in e where the pole coefficients are finite functions of the parameters ta which can be 
integrated numerically. 

Note that the coefficient functions of a certain pole, produced by the iterated decomposition of a certain 
primary sector into subsectors, are all summed over before numerical integration in order to avoid large 
numerical errors due to large cancellations. In contrast, the contributions from the primary sectors (i.e. six 
contributions in the 1 — s- 4 case) are summed over only after the numerical integration. In this way, one 
has a strong check at hand if for topological reasons (symmetry of the considered integral) some primary 
sectors have to yield identical results, while the numerical error is not substantially increased. 

We would like to emphasize that the functions occurring in the phase space integrals are quite different 
from the functions present in any virtual multi-loop integral after Feynman parametrisation: They are not 
polynomial anymore, as they unavoidably contain square-roots (and fake singularities away from the origin, 
as discussed above). This is of course not a principal problem for the algorithm, but has to be dealt with 
appropriately. 

The 0-function and the A(^secj) term in (|5.9|l become quite involved expressions after several iterations, 
but do not present any problem for the numerical integration. 

The program consists of 3 building blocks: The first one serves to perform the iterated sector decompo- 
sition until a form of the integrand is reached where all singularities are factored out explicitly as described 
above. The corresponding algorithm has been implemented in Mathematica 46 . 

Then the subtractions and the expansion in e are carried out to obtain a set of finite functions for each 
pole coefficient. These functions are then translated into FORTRAN functions. Up to this stage, the method 
is not "numerical" yet. Only the complicated structure and the large number of these functions hampers 
their analytical integration. 

^By "composite" we mean that the denominator is a nontrivial function of parameters ti, as for example {t2 + tn + te) 
present in Rs^b^ a-nd not only a trivial factor. Note that after the transformation 15.81 a previously trivial factor like l/ig 
becomes composite. 
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The third part consists in the numerical integration of these functions, where the Monte Carlo program 
BASES is used. All FORTRAN files are generated automatically using Maple [37] and perl. Perl routines 
also serve to automate compilation, job submission and summation of the results from the different primary 
sectors, such that the only "manual" input is the matrix element and the desired expansion level in e. 



5.2 Numerical Results 

The program produces an overall number of 13 subsectors for the integral Rg, 690 subsectors for i?8,a and 
321 subsectors for Rs.b- The integration time ranges from about 10 minutes for Rq to about 24 hours for 
Rs^a if a numerical accuracy better than 1% is demanded. 

The following numerical results are obtained for the master integrals: 



i?4 

Re 

R&,a 



SriqY^^" [0.08335 + 0.81946 6 + 4.14136^ + 14.398 
Sriq^y^' [0.64498 + 7.0423e + 40.5076^ + 0{e^)] , 



38, 



SriqT'-'' 



5.0003 0.0013 65.832 151.53 



0.74986 0.00009 14.001 



e 

52.911 



37.552 + 0(6) 
- 99.031 + 0{e) 



e* e 
The reducible integral i?8,r also has been calculated numerically to complete the cross-checks: 

"0.24999 0.00025 0.82251 8.4106 



R%,r — •S'r('Z 



49.748 + 0(e) 



(5.12) 
(5.13) 

(5.14) 
(5.15) 

(5.16) 



The agreement of the numerical and analytical results is better than 1%. 



6 Summary and Conclusions 



In this paper, we studied the inclusive four-particle phase space integrals of (1 4)-parton QCD matrix 
elements. These integrals contain infrared divergences due to the soft or coUinear emission of up to two 
partons, which occur also in NNLO calculations of jet observables. We demonstrated that these phase 
space integrals can be expressed as a linear combination of only four master integrals. These four master 
integrals were computed both analytically and numerically up to their fourth order terms in dimensional 
regularisation, as required for NNLO calculations. The analytic results for them are: 



i?4 = 



St {q') 
Sr {q'r'' 



2,2-2. r5(l-e)r(2-2e) 



r(3 - 3e)r(4 - 4e) 
-l + -r + ef-12 



014) 



TT 



5^ 
6 



-91 



97r2 



45C3 



90 



Iso" 



Rg^a = Sr {(f 



R 



'8.b 



Sr (q^ 



2-2e 


5 


20^2 


126C3 






3e2 


e 


2-2e 


3 


17^2 


44C3 




4? 


12e2 


e 



0{e) 



60 



(|125) 

illis) 

(|129) 



The inclusive four-particle phase space integrals of QCD 1^4 matrix elements may in future work be 
used for the construction of infrared subtraction terms for the real contributions to NNLO jet calculations. 
In the case of two-jet production in e~^e~ collisions, the 1 — s- 4 matrix elements can serve immediately as 
subtraction terms, since the only requirements on subtraction terms are that they coincide with the full 
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matrix element in all infrared singular limits (which is clearly the case) and that they can be integrated 
analytically over the inclusive phase space (which we have demonstrated here). Constructing NNLO sub- 
traction terms for the double real emission contributions to three-jet production in e~^e~ collisions from 
1^4 matrix elements, which can serve to locally approximate the full 1^5 matrix element in certain 
limits, is a yet outstanding task. To the same aim, one could also look for a way to relate already ex- 
isting NNLO subtraction terms |16II17[[TH) to 1 ^ 4 QCD matrix elements. In this context, it should be 
pointed out that the frequently used NLO dipole subtraction terms for multi-jet production in e+e^ 
annihilation can in fact be recombined to yield 1^3 QCD matrix elements. 

In the course of the calculations presented in this paper, new analytical and numerical methods for 
handling infrared divergent phase space integrals were developed. In particular, the tripole parametrisation 
of the four-particle phase space incorporates explicitly the phase space factorisation from an (n+2)-parton 
phase space to an n-parton phase space and a factor which accounts for the singular emission, as required 
for n-jet production at NNLO. This tripole formulation has the advantage that it avoids the introduction of 
non-linear transformations of momenta, which is crucial for the analytic integration. For the construction 
of an NNLO Monte Carlo parton-level event generator, a smooth mapping of all potentially singular regions 
of phase space is essential for numerical stability. Such a mapping can for example be achieved using non- 
linear momentum transformations |17lll8j . However, the phase space parametrisations used to integrate 
the subtraction terms analytically and those used in the Monte Carlo program do not need to be the same. 

Formulating the unitarity relations between phase space and loop integrals in a form which connects 
different master integrals enabled us to deduce three of the four master integrals from known multi-loop 
integrals, such that only one of them had to be calculated explicitly. Thus we were able to use multi-loop 
techniques in order to obtain results for complicated phase space integrals. 

The iterated sector decomposition, a method originally developed for infrared divergent loop integrals, 
where the poles are isolated automatically and their coefficients are calculated numerically, was extended 
here to phase space integrals. This technique is particularly powerful for (1 — > n)-particle reactions, where 
the phase space integrals depend only on one overall scale. As a consequence, the pole coefficients are 
just numbers which can be calculated numerically once and for all. This method can also be useful in 
cases where an analytic calculation of the phase space integrals cannot be achieved anymore, not only for 
inclusive cross sections, but also in view of the construction of an NNLO Monte Carlo program, where 
experimental cuts acting on the phase space boundaries can complicate the analytic integrations. 

The excellent agreement between analytical and numerical results for the four master integrals evaluated 
here provides a strong check on both calculational approaches. 
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A Massless Multi-parton Phase Space 

The n-particle phase space in the (1 —> ri)-particle decay reaction kinematics {q ^ pi + . . . + Pn) reads in 
dimensional regularisation with d = 4 — 2e space-time dimensions 



rj^{2nrS%-p,-...-Pn). (A.l) 



In the present work, we are in particular concerned about the two-, three- and four-parton phase space. 

In practical applications, it turns out that the above expression for the multi-particle phase space 
in terms of momenta is not convenient. We therefore use expressions in terms of kinematic invariants 
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Sij — 2pi ■ pj and angular volumes instead. Further, we factorise out overall rotations of the coordinate 
frame in d dimensions into the d-dimensional hypersphere dfi^ with volume 

The two-particle phase space in these variables simply reads 

dPS2 = (2^)2-'^(si2)^ dsi2 S{q^ - S12) . (A.3) 

For the three-particle phase space, one finds 

dPSs = {27rf-^'^ 2-^-"^ (g^)^ dOd-ldl7d-2 (S12S13S23)'^ dsi2 dsi3 dS23 S {q^ - S12 - Sl3 - S23) . 

(A.4) 

Finally, the four-particle phase space is 

dF54 = (27r)4-3'i(q2)3-i2l-2'i(-A4)^ e(-A4) 5{q^ - ,Si2 - S13 - S14 - .S23 - ^24 " S34) 

AVld-1 drJd-2 driti_3dsi2 dsi3dsi4ds23ds24ds34, (A. 5) 
where the Gram determinant A4 is given by 

-I-A4 = A(si2S34,si3S24, S14S23) , A(a;,y, z) + + - 2(a;y -f -I- yz) . (A. 6) 
In several parts of this paper, we also need the total n-particle phase space volume 



Pn= j dPSr,. (A.7) 

Using the unitarity relation of Section 14.31 we find a compact expression for arbitrary n: 



_ r,5-4n-2£+2nis ^3-2n-e-|-ne {\ e) / 2\n-2-|-e-ne /a 

r((n-i)(i-e))r(n(i-e))^'^^ ' ^ ■' 



in particular: 



P2 = 2-3+2^ -"'■'^^^^('^'r, (A.9) 

P3 = 2-^+^^ ^-3+2^ , ^'^^"/^ 7(9')'"'% (A.IO) 

^ r(2-2e)r(3-3e) ^ ' > 

^ ^-11+6. -5+3e ^^^^'/^ (A.ll) 

^ r(3 - 3e)r(4 - 4e) ^ ^ ' 
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